Statement of the problem

Our goal is to build a framework to study some discrete nonlinear dynamical systems and their chaotic behaviour from a numerical point of view using the Julia Programming Language. We will

1. apply this framework to the logistic map: $x_{n+1}=4r(1-x_n)x_n$ as a test
2. study two alternative unimodal maps: the tent map, and the sine map. We will observe similarities and differences
3. verify which of these characteristics are preserved in the case of a non-unimodal map: a bimodal map which we will construct as two side-by-side logistic maps.

Notes:

1. The visualizations of the bifurcation diagrams are shown in a lightweight version, where non-periodical orbits are not included. To achieve more visually appealing and resource-intensive representations, set the include_non_periodic keyword argument to true, or delete it altogether, since its default value is true.

2. The function to Lyapunov_average implements different ways of computing the Lyapunov exponent, accessible through the mode keyword argument: the "numerical" mode is based on the formula: $$ \lambda = \lim_{n\to\infty}\frac 1n\sum_{i=0}^{n-1}\ln\left|\frac{\Delta x_{i+1}}{\Delta x_i}\right|, $$ while the "analytical" mode relies on the analytical calculation of the derivative of the map, through the formula: $$ \lambda = \lim_{n\to\infty}\left(\frac{1}{n}\sum_{i=0}^{n-1}\ln|f'(x_i)|\right). $$ These limits are true only for convergence to periodic orbits, but they can be extended to chaotic orbits with some limitations.

   While the "numerical" mode is more flexible, it has the disadvantage that the computation of the Lyapunov exponent has to be optimized differently for the case of periodic orbits (where exponential convergence is best represented on the long run) vs chaotic ones (where exponential divergence is only realistic in the first steps due to fact that the orbit itself is limited).

   "analytical" mode, on the other hand, has the disadvantage that each map has to be accompanied by its analytical derivative, but the advantage that, since derivatives are oblivious of the orbit's limitedness, the implementation is stable on the long run, also for chaotic orbits. The basic implementation I propose is unable to handle the difference for the numerical case in an efficient way, so I chose to use the analytical method throughout this notebook.

Basic functions

The functions to compute relevant quantities regarding the dynamical systems are defined in a separate file to keep this one more readable. They are imported through the Julia module system and their documentation is still accessible:

(@isdefined LogisticMaps) || include("logistic_maps.jl") # Only import the module once, otherwise there will be conflicts
using .LogisticMaps

?LogisticMaps.maps #The dictionary is not exported, so we must reference it with the module name

Dictionary containing all the available maps.

Implemented keys: "logistic", "tent", "sine", "bimodal".

?logistic_map #The function is exported, so its name is available for direct use

search: logistic_map

logistic_map(x₀, r, N_steps; map="logistic")

Run the dynamical system associated to a map from maps (defaults to "logistic"), starting from point x₀, with parameter r, for N_steps discrete steps. Return a vector containing all visited positions.

Test Case: Logistic Map

Convergence at equilibrium

We can analytically compute that, for $r=1/3$ and $r=1/2$, $1/4$ and $1/2$ will respectively be equilibrium points. We can see numerically that they are stable:

using Random
using Plots

r₁= 1/3
x₀₁= rand(MersenneTwister(1233), 10) .* .74 .+ 1/4

r₂ = .5
x₀₂= rand(MersenneTwister(1233),10) .* .49 .+ 1/2;

p = plot(legend=false)
map(x-> plot!(p,x), vcat(logistic_map.(x₀₁, r₁, 15), logistic_map.(x₀₂, r₂, 15)))
plot!(1:15, ones(15) * 1/4)
plot!(1:15, ones(15) * 1/2)

display(p)

@gif for x in x₀₁ 
    draw_cobweb(logistic_map(x, r₁, 15),r₁)
end fps = 5

[ Info: Saved animation to d:\Università\Magistrale\Computazionale\esame\tmp.gif

Fixed points and periodic orbits

We can study the presence of fixed equilibrium points or periodic orbits and plot them as a function of the parameter $r$:

points = vcat([fixed_points_as_tuples(r, n=1000, map="logistic", include_non_periodic=false) for r in 0.2:.0005:1]...);

using Plots
using InspectDR
plt = scatter(points, ms=.2, ma=.2, title="Periodic orbits of the logistic map")

We can see that the stable equilibrium point, initially at $x=0$ (I'm not showing that part for visualization reasons. It's clearly visible, e.g. in ), shifts to higher values ($x=0$ remains an unstable equilibrium point), until it starts following a bifurcation pattern that quickly leads to chaos, followed by isolated windows of stable equilibrium which quickly undergo bifurcation with the same pattern.

The succession of multiplicities of stable equilibrium points looks highly unpredictable: stability regions have highly irregular extensions, so there is no hope of thoroughly characterizing their succession with simple numerical methods. Since bifurcation is poorly described by the actual measured periodicity of the system (see following plot: numerical errors lead to single points being counted twice), we can resort to the Lyapunov exponent as a tool to study this situation.

using Plots
plot(0.1:0.001:1, find_period.(0.1:0.001:1) .|> x-> x[1], legend=false, ylabel = "period", xlabel="\$r\$", title = "Periodicity of the system")

Lyapunov exponents and bifurcation points

To characterize the chaos of the dynamical system, we can compute the Lyapunov exponent as a function of the parameter $r$:

δr = .0002
r_range = 0.2:δr:1
exponents = map(r -> Lyapunov_average(0.1:.001:.99, r, start=10, map="logistic"), r_range);

scatter(r_range, exponents, ms=.5)

As the Lyapunov exponent approaches $0$, the stability of the equilibrium points or periodic orbits weakens. This suggests that the presence of peaks in the negative regions of Lyapunov exponents correspond to bifurcation points. We can detect these peaks by numerical methods, and thus calculate the position of the bifurcation points with some desired precision $\delta r$. We wrap everything up into a function so that we can use it also for other maps:

function findnegativelocalmaxima(array; ε=0.)
    result=[]
    for i in eachindex(array)[begin+1:end-1]
        if array[i-1] ≤ array[i] ≥ array[i+1] && array[i] ≤ ε
            push!(result, i)
        end
    end
    result
end

"Study the bifurcation points of the dynamical system, and plot them against"
function study_bifurcations(exponents, r_range, plt; ε=0)
    #Copy the fixed points' plot
    plt2=plot(plt)

    bifurcation_points = r_range[findnegativelocalmaxima(exponents, ε=ε)]
    for r in bifurcation_points
        plot!(plt2, [(r,0),(r,1)], label=false)
    end
    display(plt2)
    return bifurcation_points
end

study_bifurcations

bifurcation_points = study_bifurcations(exponents, r_range, plt)

5-element Vector{Float64}:
 0.2502
 0.749
 0.862
 0.886
 0.8912

This method works well at least for the first four bifurcation points. Note that numerical noise in the Lyapunov exponents may translate into artifacts in the calculations of the local maxima due to spurious oscillations. Therefore, an increase in the precision of the bifurcation point's localization (smaller $\delta r$) must correspond to an increase in the number of paths on which the Lyapunov exponent is averaged. I found that a value of $\delta r= 2\cdot 10^{-4}$, corresponding to three significant digits on the bifurcation points, gives the optimal compromise between speed and precision. Smaller values (e.g. half) cause artifacts or require more paths for the Lyapunov average. Furthermore, the calculated values for the exponents significantly depend on the start and stop parameters, so we must take our results with a grain of salt.

In any case, I wasn't able to calculate the position of the fifth bifurcation point. This makes sense since the distance between bifurcation points decreases roughly exponentially, as we can start seeing already from the first points:

# Compute the distances between bifurcation points
lengths = diff(bifurcation_points)

# Semilog plot of the distances
scatter(lengths, yscale=:log10, ylabel="distances", label = "distances")

# Find the exponent by linear fit on the logarithm
using EasyFit
fit = fitlinear(1:4,log.(lengths))

------------------- Linear Fit -------------

Equation: y = ax + b

With: a = -1.52399737530223
      b = 0.8438145313869172

Pearson correlation coefficient, R = 0.999963256887351
Average square residue = 0.00021335770320741156

Predicted Y: ypred = [-0.6801828439153128, -2.204180219217543, ...]
residues = [0.01536722126094281, -0.02381275894774637, ...]

--------------------------------------------

display(plot!([(1, exp(fit.a+fit.b)),(4,exp(4fit.a+fit.b))], label = "exponential fit"))

println("The exponent is: $(exp(-fit.a))")

The exponent is: 4.590538673751068

Other unimodal maps

We study two additional unimodal maps: the sine map, $$ x_n = r \sin(\pi x) $$ and the tent map, $$ x_n = \begin{cases}2r x & 0\leq x \leq 0.5\\ 2r(1 - x) & 0.5<x\leq 1.\end{cases} $$ Note that we redefined the classical tent map with a factor of $2$, so that for all maps the range for $r$ is $(0,1)$. Let us carry on the same study as for the logistic map, and spot any similarities or differences between the three cases.

points_sine = vcat([fixed_points_as_tuples(r, n=1000, map="sine", include_non_periodic=false) for r in 0.7:.0002:1]...);
points_tent = vcat([fixed_points_as_tuples(r, n=1000, map="tent", include_non_periodic=false) for r in 0.1:.005:1]...);

using Plots
using InspectDR
display(plt)
plt_sine = scatter(points_sine, ms=.2, ma=.2, title="Periodic orbits of the sine map", xlabel="\$r\$", ylabel = "\$x\$", legend=false)
plt_tent = scatter(points_tent, ms=1, title="Periodic orbits of the tent map", xlabel="\$r\$", ylabel = "\$x\$", legend=false)

display(plt_sine)
display(plt_tent)

As for similarities, we immediately see that all three maps have a region where $x=0$ is a stable equilibrium point. This equilibrium then shifts toward a higher value before bifurcation occurs.

In the case of the sine map, there is an incredible resemblance to the sine map. The bifurcation diagram looks like a deformed version of the logistic map's diagram. The same regions of stability and chaos are present and the bifurcation pattern appears identical.

In the case of the tent map, though, there seems to be no bifurcation at all, and, after a discontinuous shift of the equilibrium point away from $x=0$, no equilibrium is found at all.

These observations can be seen more clearly through the Lyapunov exponents:

exponents_sine = map(r -> Lyapunov_average(0.1:.001:.99, r, start=10, map="sine"), r_range);
exponents_tent = map(r -> Lyapunov_average(0.1:.01:.99, r, map="tent"), 0.1:.005:1);

using Plots
display(scatter(r_range, exponents_sine, ms=.5, xlabel="\$r\$", ylabel = "\$\\lambda\$", legend=false))
display(scatter(0.1:.005:1, exponents_tent, ms=.5, xlabel="\$r\$", ylabel = "\$\\lambda\$", legend=false))

The tent map shows no local maxima for the Lyapunov exponent in the negative region, indicating the absence of bifurcation. On the other hand, the sine map shows again the same pattern as the logistic map. By studying the bifurcation points, we can see that, though their position and difference varies from the logistic map, their distance decreases again exponentially, and with roughly the same exponent: it is like the bifurcation diagram was stretched while maintaining its proportions. This time, also the fifth bifurcation is found, and the exponential separation can be better viewed and quantified:

bifurcation_sine = study_bifurcations(exponents_sine, r_range, plt_sine)

6-element Vector{Float64}:
 0.3162
 0.719
 0.833
 0.8588
 0.8642
 0.8652

plot!(xlims=(.855,.87), ylims=(.34,.4))

lengths_sine = diff(bifurcation_sine)

# Semilog plot of the distances
scatter(lengths_sine, yscale=:log10, ylabel="distances", label = "distances")

# Find the exponent by linear fit on the logarithm
using EasyFit
fit_sine = fitlinear(1:5,log.(lengths_sine))

------------------- Linear Fit -------------

Equation: y = ax + b

With: a = -1.504667982840361
      b = 0.7405310808620347

Pearson correlation coefficient, R = 0.998614926872471
Average square residue = 0.012569473249220331

Predicted Y: ypred = [-0.7641369019783264, -2.2688048848186875, ...]
residues = [0.1451782161594034, -0.09724805423104588, ...]

--------------------------------------------

display(plot!([(1, exp(fit_sine.a+fit_sine.b)),(5,exp(5fit_sine.a+fit_sine.b))], label = "exponential fit"))

println("The exponent is: $(exp(-fit_sine.a))")

The exponent is: 4.502658422224681

Note: for this bifurcation process, we can construct the succession of the ratios between neighouring bifurcation points' distances. It has been proven that this succession has the same limit for all unimodal maps, called the first Feigenbaum constant: $$ \delta = 4.669201\dots $$ In other words, the points in the previous plot asyntotically tend to align on a straight line whose pendence is the logarithm of the first Feigenbaum constant, and this holds for all maps. Since the tent map shows no bifurcation, it can somehow be regarded as a degenerate case.

Bimodal map

I will now investigate the case of a bimodal map, i.e. a map which has two maxima. The map I chose is given by the law: $$ x_{n+1} = \begin{cases}4r (1-2x_n)\cdot 2x_n & x_n\leq 1/2\\ 4r(1-(2x_n-1))\cdot(2x_n-1) & x_n > 1/2.\end{cases} $$ Intuitively, it is a logistic map repeated twice.

plot(LogisticMaps.maps["bimodal"].(range(0,1,1000),.2), legend=false)

The first thing to look at are periodic orbits:

points_bimodal = vcat([fixed_points_as_tuples(r, n=1000, map="bimodal", include_non_periodic=false) for r in 0.01:.001:1]...);

using Plots
using InspectDR
plt_bimodal = scatter(points_bimodal, ms=.2, ma=.2, title="Periodic orbits of the bimodal map", xlabel="\$r\$", ylabel = "\$x\$", legend=false)

Which is a somewhat unexpected result: the bimodal map behaves exactly the same way as the unimodal logistic map for values of $r\in(0,1/2)$, which makes sense since in that interval one can map $x\mapsto y=2x$ and $r\mapsto s=2r$ and obtain the same form as the logistic map with half the domain both in $x$ and $r$. The second part, though, shows an initial chaotic region, with some small regularity windows, followed by a strong return to a single point stability, with bifurcations occurring at a rather slow pace. Afterwards, instability rules again.

I shall compare the two stability regists. Since it is analytically clear that the first region is identical to the logistic map case, I shall only study the second region. I will focus on the Feigenbaum constant, and test whether it is the same as the one for unimodal maps.

bimodal_range = 0.8:.0001:.86
exponents_bimodal = map(r -> Lyapunov_average(0.1:.00005:.99, r, map="bimodal"), bimodal_range);

using Plots
display(scatter(bimodal_range, exponents_bimodal, ms=.5))

bifurcation_bimodal = study_bifurcations(exponents_bimodal, bimodal_range, plt_bimodal, ε=.02)

4-element Vector{Float64}:
 0.8078
 0.8427
 0.8505
 0.8523

lengths_bimodal = diff(bifurcation_bimodal)

# Semilog plot of the distances
scatter(lengths_bimodal, yscale=:log10, ylabel="distances", label = "distances")

# Find the exponent by linear fit on the logarithm
using EasyFit
fit_bimodal = fitlinear(1:3,log.(lengths_bimodal))

------------------- Linear Fit -------------

Equation: y = ax + b

With: a = -1.4823500821532967
      b = -1.8782560387401865

Pearson correlation coefficient, R = 0.9999805517446301
Average square residue = 5.69814659682636e-5

Predicted Y: ypred = [-3.3606061208934834, -4.84295620304678, ...]
residues = [-0.005337671119729048, 0.01067534223980715, ...]

--------------------------------------------

display(plot!([(1, exp(fit_bimodal.a+fit_bimodal.b)),(3,exp(3fit_bimodal.a+fit_bimodal.b))], label = "exponential fit"))

println("The exponent is: $(exp(-fit_bimodal.a))")

The exponent is: 4.403281604541697

which, as far as we can tell, could be in good accordance with the Feigenbaum constant for unimodal maps. This suggests that there could be a transformation, affine in $r$, that could recover the form of the logistic map. Unfortunately, I was unable to find such a transformation analytically.

Additional images

Logistic map

Tent map

display.(
    draw_cobweb.(
        map(
            r -> logistic_map(.3, r,200, map="tent"), 
            [.3,.6,.9]
            ),
        [.3,.6,.9],
        start=1, map="tent"
        )
    )

3-element Vector{Nothing}:
 nothing
 nothing
 nothing